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An algorithm to create difficult Sudoku puzzles is proposed. An Ising spin-glass like Hamiltonian 
describing difficulty of puzzles is defined, and difficult puzzles are created by minimizing the energy 
of the Hamiltonian. We adopt the replica exchange Monte Carlo method with simultaneous tem- 
perature adjustments to search lower energy states efficiently, and we succeed in creating a puzzle 
which is the world hardest ever created in our definition, to our best knowledge. (Added on Mar. 
11, the created puzzle can be solved easily by hand. Our definition of the difficulty is inappropriate.) 



I. INTRODUCTION 



Sudoku, which is also called Number Place, is a kind 
of pencil puzzles [l[ . Each Sudoku puzzle has 9x9 grid, 
which has nine 3x3 subgrids. There are some numbers 
in cells. The objective of the puzzle is to complete the 
grid by filling numbers in empty cells so that each row, 
column, and subgrid contains all of numbers from 1 to 9. 
Sudoku puzzles are now extremely popular in the world, 
and recently Sudoku have attracted much attention as 
mathematical and physical points of view. In 2002, solv- 
ing Sudoku puzzle is proved to be NP-complete prob- 
lems 0. All possible Sudoku solutions is enumerated 
by Felgenhauer and Jarvis Solving Sudoku puzzles 
corresponds to find the grand state of the antiferromag- 
netic 9-state Potts model with special interactions. The 
similarity between the Sudoku problems and spin-glass 
systems has been pointed out |[. Williams and Ack- 
land defined a Sudoku Hamiltonian, and observed ther- 
modynamic phase transitions by utilizing Monte Carlo 
(MC) simulations Q. They also pointed out that the 
energy landscape of the Sudoku Hamiltonian is rugged, 
and the model show similar behavior to spin-glass sys- 
tems. From the view point of the computing, it is rather 
easy to find the solution of the given puzzle. However, it 
is not trivial to make difficult puzzles by using computers. 
In 2010, Dr. Arto Inkara created a difficult Sudoku puz- 
zle (Inkara2010) which is shown in Fig. [T](a). Later, 
he created a more difficult one in 2012 (Inkara2012) [7] 
which is shown in Fig. [T] (b). To our best knowledge, 
Inkara2012 is the world hardest Sudoku puzzle ever cre- 
ated. The purpose of the present manuscript is to create 
a Sudoku puzzle which is more difficult than Inkara2012, 
and consequently is the world hardest, by utilizing a MC 
method. 



This manuscript is organized as follows. A method 
to solve Sudoku puzzles and a definition of difficulty are 
described in Sec. [TTJ The algorithms for finding difficult 
puzzles are explained in Sec. IIIII The numerical results 
are given in Sec. IIVI and a summary and a discussion of 
further issues are given in Sec. [Vj 
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FIG. 1: Sudoku puzzles created by Dr. Inkara. (a) The puzzle 
made in 2010 (Inkara2010), which depth is 5, normal width 
is 173, and average width 179 ± 3.25, respectively, (b) The 
puzzle made in 2012 (Inkara2012), which depth is 8, normal 
width is 3599, and average width = 2257 ± 25.7, respectively. 



II. DEFINITION OF DIFFICULTY 

Difficulty of Sudoku puzzles depends on a method to 
solve them. Therefore, we describe the method to solve 
Sudoku puzzles. While there are many kinds of tech- 
niques, we adopt only two of them, pencil marks and 
recursive backtracking (sec Fig. [2]). 

Pencil Marks: Pick up an empty cell. Check all the 
numbers in the row, column, and subblock to which the 
cell belongs. Then list up all numbers which are still 
possible in the empty cell. These numbers are called 
pencil marks. If a cell has only one pencil mark, then the 
mark is the value of the cell. Repeat the above procedure 
until all empty cells have two or more pencil marks. 

Recursive Backtracking: Pick up the cell which has 
the smallest number of pencil marks. Then choose one 
of the pencil marks and assume that it is the value of the 
cell, and continue to solve the problem recursively. If the 
assumed value does not lead to the complete solution, 
then choose another value of the pencil marks. 

The algorithm to solve Sudoku puzzles involves the it- 
eration of pencil marks and recursive backtracking. Since 
the computational costs of pencil marks are cheap, we do 
not includes pencil marks to difficulty. We define three 
kinds of properties for Sudoku puzzles, depth, normal 
width, and average width, respectively. The depth is the 
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FIG. 2: (Color online) Pencil marks in sudoku cells. Each 
empty cell contains small numbers which are possible inside 
the cell. The small numbers are called pencil marks. If an 
empty cell has only one pencil mark, that mark is the value of 
the cell. In this puzzle, the smallest number of pencil marks in 
a cell is two. There are four cells which have two candidates. 



smallest number of recursive backtracking to solve a puz- 
zle. If a puzzle is solved only with pencil marks, then the 
depth of the puzzle is zero. The width is the number of 
recursive backtracking required to confirm that the puz- 
zle has a unique solution. It is possible that there are 
two or more empty cells which have the same number of 
pencil marks. In Fig. [2l the smallest number of pencil 
marks in cells is two, and there are four cells which have 
two pencil marks. A value of width depends on a choice 
of cells to perform recursive backtracking. We define two 
kinds of width, the normal width and the average width. 
Normal width: If there are two or more candidates of 
cells, chose a cell in top-left order. The choice of cells 
is deterministic, and the value of width is determined 
uniquely. Average width: If there are two or more 
candidates of cells, chose a cell randomly. The value of 
the average width is stochastic. 

The recursive backtracking process constructs a tree 
structure. The root node of the tree is the given puzzle. 
When a empty cell is chosen for recursive backtracking, 
then the number of pencil marks in the cell corresponds 
to the number of edges of the nodes. Each child node 
corresponds to a grid with an assumed number in the 
chosen cell. The tree is constructed recursively for chil- 
dren. Since we have a choice of cells to perform recur- 
sive backtracking, there are many possibility for the tree 
structure. The depth is the shortest path from the root 
node to the node describing answer between all possi- 
ble trees. The width is a number of nodes in the tree 



graph. Computational cost to search the solution is pro- 
portional to width of the graph. Therefore, we adopt 
width as difficulty of Sudoku puzzles. Since the value of 
the normal width depends on orientation of a puzzle, the 
average width is more appropriate for difficulty of Sudoku 
puzzles than the normal width. Since a number of possi- 
ble trees increases exponentially as a number of recursive 
backtracking increases, it is impossible to enumerate all 
possible trees. Therefore, we estimate the average width 
by MC sampling. In the present manuscript, we calcu- 
late average width for each puzzle from 100 independent 
samples. 

III. METHOD 
A. Model 

Suppose a solution of a puzzle, i.e., all of cells arc filled 
with numbers, is given. To create a puzzle, we have to 
remove some of numbers from the grid. The rest numbers 
are hints of this puzzle. We label cells from 1 to 81, and 
describe the state of i-th cell by a spin s^; the i-th cell 
is empty when s; = and the cell keeps the number 
of the answer when Si = 1 . A set of spin configuration 
{sj} denotes a puzzle. We define an Ising spin-glass like 
Hamiltonian of this system as 

81 

H{{ Si }) = -JU{{ Si }) + hJ2 s h (1) 

i=i 

where U denotes the internal energy given by the config- 
uration {si}, J(> 0) is the interaction energy, and h(> 0) 
is the amplitude of the external field, respectively. The 
interaction energy U is defined so that the energy of the 
system decreases when the puzzle is more difficult. Wc 
define two kinds of internal energy, depth energy Ud and 
width energy U w . Consider a puzzle which depth is d and 
normal width is w. Then the depth and width energies 
are defined to be 

U d = d, (2) 
U w = \og{w). (3) 

Since J > 0, the energies decrease as depth or width in- 
crease. We adopt logarithm for the width energy since 
width increases exponentially as a number of recursive 
backtracking increases. The energies depend on spin con- 
figuration {si}, but the relation between them is highly 
complicated. While one can easily calculate the energies 
from the given spin configuration, it is almost impossible 
to find the grand state, i.e., to find the spin configuration 
which describes the most difficult puzzle for the given so- 
lution of the puzzle This property is similar to that 
of spin-glass models. 

The second term in the right-hand side of Eq. (fTJ) in- 
crease energy when a number of hints of a puzzle in- 
creases. This term corresponds to an external field. 
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While we have to minimize only the internal energy for 
the purpose to search a difficult puzzle, we added this 
term since to apply a bias to reduce a number of hints. 
Since the spins are all-up in the initial states, it is difficult 
to decrease energy without this term. Our purpose is to 
create a puzzle width a large value of average width. But 
we find that it is difficult to obtain large values of width 
from the initial state where all spins are up. Therefore, 
we first perform MC simulations in the depth-first order, 
then we switch to that in the width-first order. 

In order to update a configuration of spins, we adopt 
the Markov Chain Monte Carlo method. Choose a spin 
randomly and flip it with the Metropolis criterion with 
the Boltzmann weight @. Note that, a spin-flip from 
down to up always increases energy, and that from up to 
down always decreases energy. If a new configuration has 
two or more solutions, then the trial is rejected. 



B. Adjusting Temperature Set 

Since the model is similar to Ising spin-glass models, 
there are many local minima in the energy land scape. In 
order to search lower-energy states efficiently, we adopt 
replica-exchange Monte Carlo (REMC) method fulfill], 
which is also called Parallel Tempering method. The 
REMC method is proposed by Hukushima and Nemoto 
in order to study hardly-relaxing systems such as spin- 
glass systems efficiently. In the REMC method, many 
replicas sharing the identical Hamiltonian are simulated 
simultaneously and temperatures of replicas are some- 
times exchanged. The REMC method requires a tuned 
set of temperatures to work efficiently. The set should 
includes a temperature which is high enough to escape 
from any local minima and a temperature which is low 
enough to search the grand state. Additionally, a num- 
ber of temperatures should be sufficient so that exchange 
ratios between adjacent temperatures are high enough. 

Usually, a temperature set is determined by prelimi- 
nary simulations and the set is fixed throughout simula- 
tions. However, many temperatures, and consequently, 
many replicas are required for Sudoku problems since 
the range of energy is wide. Therefore, we adjust the 
temperature set simultaneously throughout simulations 
to keep exchange ratios between replicas. While we can- 
not obtain the canonical ensemble without a fixed set of 
temperatures, it is not problem since we are interested 
only in the configuration having the lowest energy, 

In order to obtain a temperature set which achieves 
same exchange ratio between neighboring temperatures, 
the following procedure is proposed [Tol Il2| . 
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where M is a number of replicas, is the inverse tem- 
perature of the TO-th replica at the n-th exchange of tem- 
perature, and pj^ is the acceptance ratio of exchange be- 
tween the m-th and the (m + l)-th replicas, respectively. 
The value of p m is estimated from MC sampling between 
exchange. After convergence, the acceptance ratios will 
share the identical value as pi = p% = • ■ • = Pm—i- How- 
ever, with a large number of iterations, the temperature 
set can converge into the following trivial state 
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(7) 
(8) 



In order to prevent the trivial convergence, we determine 
a desired value of acceptance ratio p as follows, 



gn+l = o 
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(9) 
(10) 



Let N s is a number of MC steps between exchange pro- 
cesses, i.e., a number of samples to estimate the accep- 
tance ratios {p„}- If energy difference between adjacent 
replicas is extremely large, then the acceptance ratio be- 
tween the replicas becomes extremely small. Then the 
acceptance ratio p^ is estimated to be zero. Once it 
happens, we have B m = f3 m +i throughout the simula- 
tions which means that the number of replicas virtually 
decreases. In order to avoid the above, we adopts 1/N S 
instead of pj^ if the exchange is not performed between 
m-th and (m + l)-th replica in N s MC steps. Finally, we 
obtain the following procedure to adjust temperatures as, 
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The above procedure guarantees that the all tempera- 
tures have different value and a temperature increases 
when some of replicas are trapped in a local minimum. 
We choose the initial set of temperatures as 



1 



M - 1 



(14) 



The highest temperature f}\ is fixed throughout simula- 
tions. 



C. Details of Simulations 

From preliminary simulations, we adopt the interac- 
tion energy J = 100 and the external field h — 1 for 
both depth-first and width-first order calculations. First 
we create a solution of Sudoku puzzle randomly, then we 
perform the depth-first order search with the simple MC 
simulation with B = 0.05. After we find a puzzle with 
depth larger than 8, we perform the width-first search 
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FIG. 3: (Color online) Cumulative distribution functions of 
the minimum energy found by simulations with and without 
adjustments of temperatures. 2048 independent runs are in- 
vestigates for both cases. The runs with adjustments of tem- 
peratures can find lower energies than those found by runs 
without adjustments of temperatures. 



with the REMC method. We choose a number of MC 
steps between temperature exchange to be N s = 100. 
The highest temperatures is set to be /3\ = 0.01 which 
is high enough to escape any local minim. The desired 
acceptance ratio for exchange p and a number of replicas 
M are chosen to be p = 0.8 and M = 10, respectively. 
The most time consuming part of this simulation is to 
calculate energy. Therefore, we adopt the normal width 
for internal energy in Eq. ([3]) instead of the average width 
to save computational time. We list up the candidates of 
hard puzzles from simulations, and determine the hard- 
est one by calculating average width for each candidate. 
Since computational costs strongly depend on tempera- 
ture, it is inefficient to perform parallel computation for 
REMC due to load imbalance. Therefore, we assign all 
replicas to one process, i.e., a simulation of each replica 
is performed serially. 



IV. RESULTS 

In order to investigates the efficiency of the simulta- 
neous temperature adjustments, we perform simulations 
with and without the adjustments. Each simulation con- 
tains 10 replicas, and 2048 independent samples are in- 
vestigated for simulations with and without adjustments. 
We adopt the trivial parallelization, i.e., simulations are 
performed independently with different seeds of random 
numbers. Computations are carried on SGI Altix ICE 
8400EX at the Institute for Solid State Physics, the Uni- 
versity of Tokyo. Computational time is 24 hours for each 
run and the total amount of the computational time is 
about 11 CPU-core- years. 
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FIG. 4: (Color online) Typical time evolutions of a tempera- 
ture set (a) and energies of replicas (b). While ten tempera- 
tures are shown, three energies out of ten are shown for the 
visibility. The highest temperature, i.e., the lowest inverse 
temperature, is fixed to be /3i — 0.01 throughout simulations. 
One can see that the temperatures oscillate. Thanks to the 
oscillation, the replicas can escape from local minima and visit 
both high- and low-energy states. 



The cumulative distribution function (CDF) of the 
minimum energies found by the simulations are shown 
in Fig. [3] The CDF P(E) is a probability that the lowest 
energy found by each process is smaller than E. It shows 
that the simulations with temperature adjustments is 
more efficient than that without adjustments. 

Typical time evolutions of temperatures and energies 
of a run with the temperature adjustments are shown in 
Fig. SJ The reason that temperatures oscillate is as fol- 
lows. When one replica is trapped into a local minimum, 
then acceptance ratio of exchange temperature between 
the replica and its neighbor decreases. It increases tem- 
perature of the replica. After the replica escapes from 
the local minimum, then the acceptance ratio increases 
and the temperature decreases again. From the varia- 
tion rage of the lowest temperature, we estimate that 
over 35 replicas are necessary to keep the desired accep- 
tance ratio without the temperature adjustments, while 
we used only 10 replicas. It means that the simultaneous 
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FIG. 5: The hardest puzzle obtained by simulations. Its 
depth is 10, normal width is 183530, and average width is 
100571 ± 1198. 

temperature adjustments work efficiently. Figure 0] (b) 
shows that the energies of replicas fluctuate from high- 
energy states to low-energy states which means that the 
exchange of temperatures works efficiently. 

We obtain 2048 candidates of difficult puzzles from 
the simulations with the simultaneous temperature ad- 
justments. We calculate an average width for each can- 
didates, and determine the hardest one which is shown 
in Fig. [SJ The hardest puzzle's depth is 10 and aver- 
age width is 100571 ± 1198, respectively. This is about 
44 times more difficult than Inkara2012. If one adopts 
only pencil marks and recursive backtracking techniques 
to solve this puzzle, then about 50000 times recursive 
backtracking is necessary to solve it, and therefore, it is 
almost impossible to solve by hand. 

V. SUMMARY AND DISCUSSION 

We define the Hamiltonian describing the difficulty of 
Sudoku puzzles. Then creating difficult puzzles reduces 
to minimizing the energy defined by the Hamiltonian. 
We perform the REMC method to minimize the Hamil- 
tonian, and succeed to create a Sudoku puzzle which is 
much harder than Inkara2012 in our definition of diffi- 
culty. To our best knowledge, this is the world hard- 



est puzzle ever created. While the REMC is a method 
to obtain canonical ensemble of different temperatures 
simultaneously, we propose the REMC with simultane- 
ous adjustments of temperature which does not guaran- 
tee canonical ensembles. The results presented in this 
manuscript demonstrates that the REMC method is use- 
ful not only for physical problems, but for general opti- 
mization problems. 

A definition of difficulty strongly depends on solving 
algorithms. We adopt only two techniques, pencil marks 
and recursive backtracking, while there are many kinds 
of techniques to solve Sudoku puzzles. Therefore, the dif- 
ficulty defined in the present manuscript can be different 
from the actual feeling. But once the definition of the 
difficulty is given, then difficult puzzles can be created in 
that definition since the method proposed in the present 
manuscript is general and it is independent of a definition 
of difficulty. 

While we adopt the REMC method here, there are 
other optimization methods, such as genetic algorithm 
or simulated annealing, and so forth. It is one of the 
further issues to compare such algorithms with REMC. 
Since we are interested only in the lowest energy state, it 
is not necessary to achieve canonical ensemble for given 
temperature. Therefore, we have a choice of the transi- 
tion probability. While we adopt the Boltzmann weight 
both for MC and REMC, it is possible to accelerate find- 
ing lower energy states by adopting general transition 
probability [l3j . 

In the present manuscript, we create a difficult Su- 
doku puzzle. It is not always true that a difficult puzzle 
is interesting one. However, if we can define a quantity 
describing how interest a puzzle is, then we can create 
interesting puzzles by adopting the similar method pro- 
posed here. 

The programs used in the present manuscript are pub- 
lished as open source software [13] . 
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